Mon. Not. R. Astron. Soc. 000,ITHT3lO Printed 29 April 2010 (MN M£X style file v2.2) 



On the identification of merger debris in the Gaia Era 
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ABSTRACT 

We model the formation of the Galactic stellar halo via the accretion of satellite 
galaxies onto a time-dependent semi-cosmological galactic potential. Our goal is to 
characterize the substructure left by these accretion events in a close manner to what 
may be possible with the Gaia mission. We have created a synthetic Gaia Solar Neigh- 
bourhood catalogue by convolving the 6D phase-space coordinates of stellar particles 
from our disrupted satellites with the latest estimates of the Gaia measurement errors, 
and included realistic background contamination due to the Galactic disc(s) and bulge. 
We find that, even after accounting for the expected observational errors, the result- 
ing phase-space is full of substructure. We are able to successfully isolate roughly 50% 
of the different satellites contributing to the 'Solar Neighbourhood' by applying the 
Mean-Shift clustering algorithm in energy- angular momentum space. Furthermore, a 
Fourier analysis of the space of orbital frequencies allows us to obtain accurate esti- 
mates of time since accretion for approximately 30% of the recovered satellites. 

Key words: galaxies: formation - galaxies: kinematics and dynamics - methods: 
analytical - methods: TV-body simulations 



1 INTRODUCTION 

How galaxies form and evolve remains one of the most in- 
teresting and challenging puzzles in astronomy. Although 
a great deal of progress has been made in the last decade 
many questions await to be addressed. Both, theory and ob- 
servations are now converging about a key ingredient of this 
formation process; galaxies as our own Milky Way seem to 
have experienced the accretion of smaller objects that have 
come together than ks to the relentless pull of gravity (e.g. 
IWhite fc Reeslll978l ). 

From the observational side, several new studies have 
revealed the presence of a large amount of stellar streams 
in galactic stellar haloes, echoes of ancient as well as re- 
cent and ongoing accretion events. These stellar streams 
have been preferentially found in th e outer regions of galax- 
ies. The Sagittarius Tidal Streams dlbata. Gilmore fc Irwinl 
Il994l : llbata et al.ll20oibl:lMaiewski et al.ll2003h and the Qr- 
phan Stream ( Belokurov et alJ 20071 ) are jus t two examples 



|200S| . I2009T ) are an unambiguous proof that accretion is in- 
herent to the process of galaxy formation. 

On the theoretical side, models of the formation of stel- 
lar haloes in the ACDM cosmogony hav e been able to explain 
their gross structural properties (e.g. iBullock fc Johnston! 
120051 : iDe Lucia fc Helmill2008l : ICooper et alj|20ld) . In these 
models, the inner regions of the haloes (including the So- 
lar Neighbourhood) typically formed first and hence contain 
information about the most ancie nt accretion events that 
the galaxy has experienced (e .g. White & Springcl 2000l : 
Helmi. White fc Springell 120031 : bullock fc Johnstonl 120051 : 



Tumlinsonll2010l) . Previous studies have predicted that 



of satellite debris in the Milky Way (see al s olNewberg et al 
2002; llbata et al]|2003l: lYannv et all 120031: iBelokurov et al 



20061 : iGrillmain 120061 : IStarkenburg et al.ll2009l h The abun- 



dant substructures found not only in the hal o of M31 (e.g. 
llbata et al.l l2001al : iMcConnachie et al.l [2009T). but also in 
several haloes of other galaxies (e.g. Martfnez-Delgado et al.l 
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eral hundreds kinematically c old stellar streams sh ould be 
present in the Solar Vicinity l|Helmi fc White! 1 19991 ). How- 
ever, the identification of these streams is quite challeng- 
ing especially because of the small size of the samples 
of halo stars with accurate 3-D velocities currently avail- 
able. Full phase-space information is necessary because of 
the very short mixing time-scales in the inner regions of 
the halo. Some progress has been made towards build- 
ing such a catalog ue but only a few detections have been 
reported to date (|Helmi et all 1 19991 : iKlement et~aH 120091 : 
ISmith et al1l2009h. These have made use of c atalogues such 
as SE GUE l|Yannv et all I2009T ) and RAVE dZwitter et all 
120081) in combination with Tyc ho (|H0g et alJl200C| ) and Hir> 
parcos (|Perrvman et al.lll997l h Clearly, this field will only 
advance significantly with the launch of the astrometric 
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satellite Gaia (|Perrvman et al.ll200ll ) . This mission from the 
European Space Agency will provide accurate measurements 
of the 6-D phase-space coordinates of an extraordinarily 
large number of starfl. Together with positions, proper mo- 
tions and parallaxes of all stars brighter than V = 20, Gaia 
will also measure radial velocities down to a magnitude of 
V = 17. 

To unravel the accretion history of the Milky Way re- 
quires the development of theoretical tools that will ulti- 
mately allow the identification and characterization of the 
substructure present in the Gaia data set. Several authors 
have studied the suitability of var ious spaces to isolate de- 
bris from accretion events (see, e.g. Hclmi & de Zceuw 2000; 
Knebe et al.ll2005l;]Arifvanto fc Fuchsl2006l;lFont et allhoOrj ; 



Helmi et al l 120061: ISh arma fc JohnstorJl2009h. Recently 
McMillan fc Binnevl l|200Sl ), followed by iGomez fc Helmil 



(2010J, uereafter GH10), showed that orbital frequencies 
form a very suitable space in which to identify debris from 
past merger events. In frequency space individual streams 
from an accreted satellite can be easily identified, and their 
separation provides a direct measurement of its time of 
accretion. Furthermore, GH10 showed for a few idealized 
gravitational potentials that these features are preserved 
also in systems that have evolved strongly in time. While 
very promising, these studies have focused on a single ac- 
cretion event onto a host galaxy. In reality, we expect the 
Galactic stellar halo to have formed as a result of mul- 
tiple merger events (where most of its mass should have 
originated in a handful of significant contrib utors; see e.g. 
iDe Lucia fc Helmil Eiosl ; ICooper et all l2010h . As a conse- 
quence, it is likely that debris from different satellites will 
overlap in frequency space, complicating their detection. 

In this paper we follow multiple accretion events in a 
live (cosmologically motivated) galactic gravitational poten- 
tial. We focus on how the latest estimates of the measure- 
ment errors expected for the Gaia mission will affect our 
ability to recover these events. In particular we study the 
distribution of stellar streams in frequency space and the 
determination of satellite's time since accretion from stellar 
particles located in a region like the Solar Neighbourhood. 
Our paper is organized as follows. In Section 2 we present 
the details of the TV-body simulations carried out, as well as 
the steps taken to generate a Mock Gaia stellar catalogue. 
In Section 3 we characterize the distribution of debris in a 
Solar Neighbourhood-like sphere in the absence of measure- 
ment errors, while in Section 4 we focus on the analysis of 
the Mock Gaia catalogue. We discuss and summarize our 
results in Section 5. 



2 METHODS 

We model the formation of the stellar halo of the Milky Way 
using TV-body simulations of the disruption and accretion of 
luminous satellites onto a time-dependent galactic poten- 
tial. We describe firstly the characteristics of the galactic 
potential and of the population of satellites, as well as the 
TV-body simulations. Finally we outline the steps followed 
to generate a mock Gaia star catalogue. 

1 Details about the latest Gaia performance numbers may be 
found at: http://www.rssd.esa.int/gaia 



Table 1. Parameters of the present day Milky Way-like potential used 
in our TV-body simulations. Masses are in Mq and distances in kpc. 



Disc 


Bulge 


Halo 


M d = 7.5 x 10 10 


M b = 2.5 x 10 10 


M vir = 9 x 10 11 


r a = 5.4; r b = 0.3 


r c = 0.5 


r vir = 250 






c = 13.1 



In this paper we adopt a flat cosmology defined by 
£l m = 0.3 and Qa = 0.7 with a Hubble constant of H(z = 
0) = Ho = 70 km s^Mpc" 1 . 



2.1 The Galactic potential 

To model the Milky Way potential we chose a three- 
component system, inc luding a Miyamoto-Nagai disc 
jMivamoto fc Nagai|[l975l ') 



GM d 



sjR 2 + {r a + s/Z 
a Hernquist bulge (jHcrnouist 1990T ). 



2^)2 



GM b 
r + r c 



(1) 



(2) 



and NFW dark matter halo (jNavarro. Frenk fc Whitcll 19961 ) 

3-halo = -. ~ log ( 1H I . (3) 



r [log(l + c) - c/(l + c)} l0g y 1 + r t 

Table Q] summarizes the numeri cal values o f the 

parameters at redshift z = |Smith et alj 120071 ; 

ISofue. Honma fc Omodakal [2009). The circular velocity 
curve in this model takes a value of 220 km s _1 at 8 kpc 
from the galactic centre, and is shown in Figure [1] 

To model the evolution of the Milky Way potential 
we allow the characteristic parameters to vary in time. 
For the dark matter halo, the evolution of its virial mass 
and its concentration as a function of redshi ft are given by 
IWechsler et~ai1 (|2002h and lZhao etafl {2003) 



M vil (z) = M viT (z = 0)exp(-2a c z), 



(4) 



where the constant a c = 0.34 is defined as the formation 
epoch of the halo, and 



c(z) 



c(z = 0) 

1 + 2 ' 



(5) 



For the disc a nd bulge components we fo llow the pre- 
scriptions given bv lBullock fc Johnston! (120051 ). i.e., 



M d ,b(z) = M d>b (z = 0) 
for the masses and 



r a ,b,c(z) = r a , b ,c(2 = 0) 



M v i r (z) 
M vir (2 = 0) 

r v iv{z) 



r vil (z = 0) 



(6) 



(7) 



for the scalelengths. Here r vlT is the virial radius of the dark 
matter halo. Its evolution can be expressed as 



rvir (2) = 



3M vir (z) 



AnA vir (z)p c (z) 



1/3 



(8) 
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Figure 1. Circular velocity curve as a function of distance from 
the galactic centre. The thick solid lines represent the total cir- 
cular velocity used in the suite of simulations at z = (black) 
and z ~ 1.85 (grey), respectively. The individual contributions 
from the dark matter halo, the disc and the bulge at present time 
are shown by the thin dashed, dotted and dashed dotted lines, 
respectively. 



wher e A v j r (z) is the virial overdensity (|Brvan fe Normar] 

(HI), 

A vir (z) = 18tt 2 + 82[fi(z) - 1] - 39[fi(z) - l] 2 (9) 
with Q(z) the mass density of the universe, 

n m (i + 2 ) 3 



Q(z) = 



tt m (l + z) 3 + fV 



(10) 



and p c (z) is the critical density of the universe at a given 
redshift, 



Pc(z) = 



3H 2 (z[ 
8ttG 



with 



H(z) = Ho^/vTa + Q m (l + zf 



(11) 
(12) 



2.2 Satellite galaxies 

2. 2. 1 Internal properties 

We assume the properties of the progenitors of our model 
stellar halo (i.e. the satellites at high redshift) are similar to 
those at present day. The Galactic stellar halo has a total 
luminosity of ~ 10 9 Lq . To obtain a population of satellites 
that, after accretion, produces a stellar halo with this total 
luminosity, we use the luminosit y distribution functio n of 
the Milky Way satellites given bv lKoposov et all (|2008h 

^ = 10xl0°- 1(M ^ +5) . (13) 
dMv 

Figure [2] shows the luminosity function of our model satellite 
population obtained by randomly drawing 42 objects. 

For the initial st ellar structu re of the satellites, we as- 
sume a King profile |Kinel[l966l ) with a concentration pa- 
rameter c w 0.72. We derive the structural parameters for 
each of our satellites using the scaling relations 



lot 



lop 



^- - 3.53 log av _ 
Lq km s 



km s 



1.15 log 



Rc 
kpc 



2.35, 



1.64, 



as given by iGuzman. Lucev fe Bower! (1 19931 ) for the Coma 
cluster dwarf ellipticals galaxies. Here R c is the core radius 
of the King profile, and a v is the central velocity dispersion. 
We note that the total mass of the King model as defined 
by these parameters may differ from that implied by the 
satellite's luminosity L. Therefore we also implicitly allow 
the presence of dark matter in our model satellites. 



2.2.2 Orbital properties 

The density profile of the stellar halo is often parametrized 
in a principal axis coordinate system as 



p(x,y,z) = po 



v 2 z 2 

x 2 + ^- + ^r 

p z q z 



(14) 



where n is the power-law exponent, p and q the minor- and 
intermediate-to-mayor axis ratios, a the sca le radius and po 
the stellar halo density at a radius ro- (See lHelmlll2008l . for 
a complete review of recent measurements of these param- 
eters.) Here for simplicity, we shall assume a — (as the 
halo appears to be rather concentrated), p = q = 1 and 
n = —3. At the solar radius, ro = Rq = 8 kpc, po cor- 
responds to the local stellar halo density. In this paper we 
adopt po = 1.5x 10 4 Mq kpc -3 , as given bv lFuchs fc Jahreifil 
(|l998h . 

From this density profile we randomly draw positions 
at redshift z = for the guiding orbit of each one of our 
satellites. To generate their velocities we follow the m ethod 
described in section 2.1 of iHelmi fc de Zeeuwl (|200Qt) . The 
main assumption made is that the stellar halo can be mod- 
elled by a radially anisotropic phase-space distribution func- 
tion, which is only a function of energy, E, and angular mo- 
mentum, L. 

This combined set of orbital initial conditions guaran- 
tees that at z = 0, the observed stellar halo density profile 
and velocity ellipsoid at the solar radius are roughly matched 
by our model. These initial conditions are then integrated 
backwards in time for ~ 10 Gyr (z as 1.85) in the time- 
dependent potential described in Section \2. II Therefore, we 
now have obtained the set of initial positions and veloci- 
ties required for forward integration in time of our iV-body 
simulations. 

We do not include any numerical treatment of dynami- 
cal friction in our simulations. Therefore, and to mimic the 
effects of this process, we have assigned the five most bound 
orbits at present time to the five most massive satellites. For 
the rest of the satellites the orbits are assigned randomly. 



2.2.3 Numerical treatments 

The numerical simulations are, in most cases, self- 
consistently evolved using the mas sively parallel TREESPH 
code GADGET-2.0 (|Springelll2005l ). The number of particles 
used to simulate each satellite depends on its total luminos- 
ity, as explained bellow (Section I2.3|l . This number ranges 
from a minimum of 2.56 x 10 5 up to 10 7 particles. To avoid 
very large computational times we neglect self-gravity for 
the five most massive satellites. We do not expect this to 
significantly affect our results since these satellites inhabit 
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Koposov et al. 2008 
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Figure 2. Luminosity function for the model satellite galax- 
ies used in this suite of simulations (histogram). The (grey) 
dashed line shows t he 'all sky SDSS' luminosity function by 
iKoposov et~afl ll2008f) . 



the inner regions of the host (see Section r2.2.2|l and therefore 
they are rapidly disrupted. For the gravitational softening 
parameter we choose e = 0.025i? c (where R c is the core ra- 
dius of the system). After allowing each satellite to relax in 
isolation for 3 Gyr, we launch it from the apocentre of its 
orbit and let it evolve for ~ 10 Gyr under the influence of 
the time-dependent galactic potential. 



2.3 Generating a mock Gaia catalogue 
2.3.1 N-body particles to stars 

To generate a mock Gaia catalogue we need to 'trans- 
form' our TV-body particles into 'stars'. Therefore we as- 
sign absolute magnitudes, colours and spectral types to each 
particle in our simulatio ns. We use the IAC-STAR code 
l|Aparicio fc Gallart! 2004) which generates synthetic colour- 
magnitude diagrams (CMDs) for a desired stellar population 
model. As output we obtain the My and (V — I) colour for 
each stellar particle. 

We assign to each satellite a stellar population model 
l|Pietrinferni et all I2004T ) with a range of ages from 11 to 
13 Gyr, -2 < [Fe/H] < - 1.5 dex, [q/FeJ = 0.3, and a 
Kroupa initial mass function (|Kroupa et al.lll993l ). All these 
parameters are consis tent with th e observed values in the 
Galactic stellar halo (Helmi 2008, and references therein). 
An example of the resulting CMD is shown in Figure [3] 

Due to the limited resolution of our TV-body simulations 
it is not possible to fully populate the CMDs of our satel- 
lites. Therefore each satellite is populated only with stars 
brighter than My « 4.5, which corresponds to an appar- 
ent magnitude V = 16 at 2 kpc from the Sun. This choice 
represents a good compromise between the numerical resolu- 
tion of our simulations and the limiting apparent magnitude 
Vii m for which Gaia will measure full 6D phase-space coor- 



Figure 3. Example of a synthetic CMD used to populate our 
model satellite galaxies with 'stars'. The dashed line indicates the 
limiting magnitude considered in our simulations (My fs 4.5). 



dinates (astrometric and photometric measurements extend 
to V = 20, but spectroscopic measurements reach V ~ 17). 

As explained in Section 12.2.11 different satellites have 
different total luminosities and, thus, the number of stars 
down to magnitude My ai 4.5 varies from object to object. 
Therefore, in each satellite a different number of TV-body 
particles is converted into stars. To est i mate this number 
we use the CMD code by iMarigo et al.l l|2008l ) , which pro- 
vides a stellar luminosity function normalized to the initial 
stellar mass of a given population. The age and metallicity 
of the stellar population model is fixed to the average values 
described above and the initial mass is set by the total lu- 
minosity of each satellite. Finally, accumulating the number 
of stars down to My « 4.5 we obtain the number of TV-body 
particles that need to be transformed into stars. This num- 
ber range from 660 for the faintest satellite (My = —5) to 
1.4 x 10 r for the brightest one (My = -15.9). 



2.3.2 Disc and Bulge contamination 

The Gaia catalogue will contain not only stars from the halo, 
but also stars associated to the different components of the 
Galaxy. These stars act as a smooth Galactic background 
which could in principle complicate the task of identifica- 
tion of debris associated to acc retion events. To take into ac- 
count this stellar background, iBrown. Velasquez fc Aguilarl 
(2005, hereafter BVA05) created a Monte Carlo model of 
the Milky Way. This model consists of three different com- 
ponents: a disc, a bulge, and a smooth stellar halo; where 
the latter was meant to take into account the set of halo 
stars formed in-situ. In this work we only consider the back- 
ground contamination by the disc and bulge components of 
the Monte Carlo Galactic model of BVA05, since we build 
our stellar halo completely from disrupted satellites. The 
disc model consists of particles distributed in space accord- 
ing to a double exponential profile. For the bulge model a 



Identifying merger debris in the Gaia era 5 




—i 

-3000 -2000 -1000 1000 2000 3000 

L [kpc km s ] 



Figure 4. Distribution of particles in E — L z space after 10 Gyr 
of evolution. Different colours represent different satellites. 

Plummer density distribution (|Plummeijll91ll ) is adopted. 
To each particle an absolute magnitude and spectral class 
is as signed according to the H ess diagram listed in table 4- 
7 of iMihalas fc Binnevl (|l98lf ). Finally, the kinematics are 
modelled assuming different velocity ellipsoids for stars in 
the disc of each spectral type whereas an isotropic velocity 
dispersion is assumed for bulge stars. For more details, see 
section 2.1 of BVA05. 

2.3.3 Adding Gaia errors 

To convolve the positions and velocities of our mock cata- 
logue of stars with the expected Gaia errors we have followed 
the steps described in section 4.1 of BVA05. Here we give 
a short overview of the procedure but we refer the reader 
to BVA05 for more details. As a first step, we transform 
Galactocentric positions and velocities to heliocentric coor- 
dinates, rhdio = r ga i — f© and Uhoiio = v ga i — Wq. Subse- 
quently, these quantities are transformed into radial veloc- 
ity and five astrometric observables that Gaia will measure, 
i.e., the Galactic coordinates (I, b), the parallax zu, and the 
proper motions (Xi» = fii cos b and fj,b- The next step con- 
sists in convolving with the expected Gaia errors according 
to the accuracy assessment described in the Gaia web pages 
at ESA (http://www.rssd.esa.int/gaia). Note that the errors 
vary most strongly with apparent magnitude, with a weaker 
dependence on colour and position on the sky. Finally, we 
transformed the particles' observed phase-space coordinates 
back to the Galactocentric reference frame for our analysis. 



3 CHARACTERIZATION OF SATELLITE 
DEBRIS 

In this section we analyse how the debris of our 42 satel- 
lites is distributed in various projections of phase-space, 




Figure 5. Distribution on the sky (Z, 6) of the stellar particles 
located inside a sphere of 2.5 kpc radius centred at 8 kpc from 
the galactic centre. Different colours represent different satellites. 



before taking into account how the Gaia observations will 
affec t this distribution. In c omparison to previous works 
fe.g. lHelmi fc de ZeeuwI feoOO), recall that our satellites have 
evolved in a live potential for 10 Gyr. 



3.1 Traditional spaces 

Figure [4] shows the distribution of ~ 2.5 x 10 4 randomly 
chosen particles from each satellite in the space of energy, 
E, and the z-component of the angular momentum, L z . The 
different colours represent different satellites. Note that in 
this projection of phase-space a large amount of substruc- 
ture is still present, despite the strong evolution of the host 
gravitational potential (e.g. the total mass increased by a 
factor ~ 3.5). 

We focus now on the stellar particles from our mock 
Gaia catalogue (prior to error convolution) located within a 
'Solar Neighbourhood' sphere of 2.5 kpc radius, centred at 
8 kpc from the galactic centre. Figures [5] and [6] show their 
distribution on the sky and in velocity space, respectively. 
Inside this sphere we find approximately 2 x 10 4 stellar parti- 
cles coming from 22 different satellites that contribute with, 
at least, 20 stars brighter than Mv ~ 4.5. The distribution 
on the sky is very smooth and, thus, disentangling merger 
debris in the 'Solar Neighbourhood' by only using spatial in- 
formation is not obvious. On the contrary, substructure can 
be observed in velocity space but this is clearly less sharply 
defined than what is found in the pseudo integrals of mo- 
tion E-Lz space (e.g. Figure [4}. We turn next to the space 
of orbital frequencies where we also expect to find much 
dumpiness. 



3.2 Frequency space 

The orbital frequencies of the stellar particles are computed 
as follows (see also section 5.1 of GH10): 

• We assume a fixed underlying potential: as described in 
Section ED at 2 = 0. 

• We integrate the orbits of each stellar particle for ap- 
proximately 100 orbital periods. 

• We obtain their orbital freque ncies by applying the 
spectr al analysis code developed by ICarpintero fc Aguilarl 
(|l998h . 
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Figure 6. Distribution in two different projections of velocity 
space of the stellar particles shown in Fig. \E\ 



3.2.1 Generalities 

The distribution of the stellar particles in frequency space is 
shown in the left panel of Figure [7] We use the same colour 
coding as in previous figures. Note that satellite galaxies do 
not appear as a single and smooth structure, but rather as 
a collection of well defined and small clumps. These small 
clumps are associated to the different stellar streams cross- 
ing our Solar Neighbourhood sphere at present time. In a 
time independent gravitational potential, streams are dis- 
tributed in frequency space along lines of constant Q r and 
Qj,. However, in a time dependent potential, such as the one 
considered here, streams are distributed along lines with a 
given curvature, depending upon the rate of growth of the 
potential (see also GH10). 

The left-hand panel of Figure [7] also shows that in the 
tt<p vs. £l r — fl^ space satellites may strongly overlap. This 
situation can be considerably improved by adding the z- 



component of the angular momentum as an extra dimension 
to this space, as shown in the right-hand panel of Figure [7] 

The distribution of debris in the space of Q r — fi</> vs. L z 
is very comparable to that in E vs. L z , as a coarse compar- 
ison between Figure U and Figure [7] will reveal. Note as well 
that as in E — L z space, from this projection it is possible, 
solely by eye inspection, to isolate a few accreted satellites. 

From this Figure we can also notice that satellites on 
low frequency orbits have a smaller amount and a set of 
better defined streams than those on orbits with high fre- 
quency. The reason for this is twofold. Firstly, satellites on 
low frequency orbits spend most of their time far from the 
centre of the host potential and therefore have longer mix- 
ing time-scales, as opposed to those on high frequency or- 
bits (short periods). Secondly, potentials as the one consid- 
ered in this work may admi t a certain degree of chaoticity 
(e.g. ISchuster fc Allenlll997l ). Satellites on highly eccentric 
short period orbits come close to the galactic centre and 
may be 'scattered' via chaos. Such chaotic orbits do not 
have stable fundamental frequencies (since they have fewer 
integrals of motion than needed and thus wander through 
phase-space). Therefore, their structure in frequency space 
is rapidly erased. 

It is important to note that our desire to analyse a 'Solar 
Neighbourhood' sphere has resulted in certain limitations. 
For example, some satellites do not contribute at all to this 
volume because of their particular orbit. In addition, most 
of the faint satellite galaxies have only a few, if any, 'stars' 
with absolute magnitude above the limiting magnitude of 
our Mock Gaia catalogue and therefore they are not 'ob- 
servable', either. 



3.2.2 Estimating the time of accretion 

An important feature of frequency space is that an estimate 
of the time since accretion of a satellite can be derived by 
measuring the separatio n between adjacent stream s along 
the Q,^ or Q, r directions (jMcMillan fc BinnevllioOol . GH10). 
This characteristic scale may be estimated through a Fourier 
analysis technique as follows 

• We create an image of the scatter plot in frequency 
space, by griding the space with a regular N x N mesh of 
bin size A and count the number of stellar particles on each 
bin. 

• We apply a 2-D discrete Fourier Transform to this image 
and obtain H{k r ,k$) 

• We compute a one-dimensional power spectrum along a 
thin slit centred on the wavenumbers k r /(NA), k</,/(NA) = 
as 

P(0) = ^\H(0,0)\ 2 , 

P(k+) = 7^ [\H{k^ 0)| 2 + \H{-kt, 0)| 2 ] 

N (15) 
for ^ = l,...,(f -1), 

P(N/2)^^\H(-N/2,0)\ 2 
and analogously for P(k r ). 



Identifying merger debris in the Gaia era 7 




60 



50 



15 20 25 

n [Gyr" 1 ] 



40 

O 



a 
i 



30 



20 



10 



























.2 


! 















-3000 -2000 



-1000 







1000 

[kpc km s -1 ] 



2000 3000 



Figure 7. Distribution of stellar particles in frequency (left-hand panel) and L z vs. Q r — (right-hand panel) space located inside a 
sphere of 2.5 kpc radius at 8 kpc from the galactic centre. The colour-coding is the same as in previous figures. Note that this distribution 
of particles was obtained without taking into account the expected Gaia observational errors. 



• We identify the wavenumber fo of the dominant peak 
in the spectra, which corresponds to a wavelength equal to 
the distance between adjacent streams in frequency space. 

• The estimate of the time since accretion is t& cc = 2nfa 
(for more details, see section 3.2.3 of GH10). 



of the stellar particles became unbouncQ. For each satellite 
we obtain a t acc of 9, 9.4 and 9.7 Gyr, respectively. In con- 
gruence with GH10, we find that this method provides a 
lower limit to the actual time since accretion, differing by 
15 — 25%. This is very encouraging since GH10 did not study 
such a complex growth of the host potential. 



We have applied this method to three different satellites 
from our simulations. Two of these satellites (labelled num- 
ber 1 and 2 in the right-hand panel of Figure [7|) have a low 
frequency guiding orbit and, just by eye inspection, can be 
isolated in L z vs. space. The third one (number 3) 

has a high frequency guiding orbit and overlaps with some 
other satellites in this space. 

In Figure [H] we zoom-into the distribution of particles 
in frequency space for each satellite separately. From this 
figure we clearly appreciate the large number of streams 
that each of these satellites is contributing to this 'Solar 
Neighbourhood' volume. Although the number of streams 
apparent in this Figure is very large (44, 30 and 59, for 
satellite 1, 2 and 3, respectively) , this is consistent with the 
models of iHelmi fc Whitel (|l999l V who predicted a total of 
~ 300 — 500 streams locally, i.e. in an infinitesimal volume 
around the Sun. This is encouraging particularly because the 
models presente d here are much more realistic than those by 
IHelmi fc Whitel ljl999h . 

To obtain an estimate of the time since accretion for our 
satellites we compute the normalized power spectrum along 
the k r — direction. The results are shown in the bottom 
panels of Figure[8] From left to right, the peak with the high- 
est amplitude in the spectrum is located at a wavenumber 
of fo = 1.24, 1.41 and 1.21 Gyr, respectively. These values 
correspond to an estimated time since accretion of 7.8, 8.9 
and 7.6 Gyr, respectively. These are in reasonable agreement 
with the true value, which we define as the time when 80% 



4 ANALYSIS OF THE MOCK GAIA 
CATALOGUE 

We will now analyse the phase-space distribution of stellar 
particles located in a volume in the 'Solar Neighbourhood', 
as may be observed in the future by Gaia. We study how 
the Gaia errors, as well as the contamination from other 
Galactic components, would affect our ability to identify and 
characterize merger debris. 

4.1 Contamination by disc and bulge 

As in Section \3. 21 we restrict our analysis to the stellar par- 
ticles located inside a sphere of 2.5 kpc radius centred at 
8 kpc from the galactic centre. To account for the maximum 
possible degree of background contamination, we consider 
all the stars for which Gaia will measure full phase-space 
coordinates (i.e., all stars brighter than V = 17), according 
to the Monte Carlo model of the disc and bulge. 

Disc particles largely outnumber those in the stellar 
halo, generally providing a prominent background. However 
much of this can be removed by considering the Metallic- 
ity Distribution Function (MDF) of the Galactic compo- 
nents. While the MDF of the halo peaks at approximately 

2 For simplicity we assume that a particle becomes unbound 
when it is found outside the initial satellite's tidal radius, ob- 
tained from the initial King profile. 
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Figure 8. The top panels show for three different satellites, the distribution of stellar particles in frequency space located inside a sphere 
of 2.5 kpc radius at 8 kpc from the galactic centre. The bottom panels show the 1-D normalized power spectra along the k r = direction, 
obtained as explained in Section lij.2.21 The wavenumber of the dominant peak (denoted by the vertical lines in these panels) is used to 
estimate the accretion time of the satellites. Other peaks in the spectra are associated to either the harmonics of this wavenumber or to 
the global shape of the particle's distribution in frequency space. 



[Fe/H] = -1.6, that of disc peaks at [Fe/H] « -0.2. As a 
consequence, stars from the disc are in general more metal- 
rich and therefore a simple cut on metallicity could be of 
great help to isolate halo stars in our sample. 

It is therefore important to characterize well the metal 
poor tail of the G alactic disc MDF. We use the model by 
llvezic et al.l ((2008) to represent this tail: 

p(x) = Gi(x|/Lti,o-i) + 0.2 G 2 (x\n2,(T2), x = [Fe/H], 

(16) 

and 

G(*j M ,a) = -^exp - (3 r/ )2 . (17) 

Here <Ti = 0.16 dex, 02 = 0.1 dex, fi2 = —1.0. llvezic et all 
(2008) found /ii to vary as a function of the height from the 
plane of the disc, but for simplicity we fix /ii = —0.71. The 
contributio n of stars mo r e met al-rich than [Fe/H] « —0.4 to 
the study of llvezic et alj (2008) is small (because of a restric- 
tion imposed on the colour range). Therefore to account for 
this population in our model, we normalize p([Fe/H]) with 
a constant a. The numerical value of a is such that at its 
peak value (located at [Fe/H] « —0.7) p([F e/H]) matches the 
MDF of the Geneva-Copenhagen survey (|Nordstrom et af] 
12004 ) . 



Having derived a MDF for our disc model, we can pro- 
ceed to apply a cut on metallicities [Fe/H] ^ —1.1 to the 
disc particles present in our 'Solar Neighbourhood' sphere 
of 2.5 kpc radius. Although in our simulation halo stars 
have metallicities between —2.0 ^ [Fe/H] ^ —1.5, stars with 
higher metallicities are known to be present in the Galactic 
stellar halo. Hence the more generous cut. This cut leads to 
a subsample of only « 9.3 x 10 3 'disc stars' brighter than 
V = 17 out of a total of 4.1 x 10 7 found in our model. 

A similar approach can be followed for the bulge com- 
ponent. In this case metallicities are assigned acco rding to 
the ob served MDF of the Galactic bulge, as given bv lZoccalil 
(2009). This MDF is modelled as a single gaussian distribu- 
tion with a mean metallicity ([Fe/H]) = —0.28 dex and a 
dispersi on a = . 4 dex , as found in the outermost field anal- 
ysed bv lZoccalil l|2009h (located at a latitude of b = -12°). 

After removing all stellar particles with [Fe/H] ^ —1.1, 
out of the 2.9 x 10 4 bulge stellar particles with V < 17 
originally present in our 'Solar Neighbourhood' sphere we 
are left with only 784 stars. This very small number suggests 
that the contamination from this galactic component should 
not strongly affect the detection of substructure in phase- 
space. 
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Figure 9. Distribution of stellar particles in frequency (left-hand panel) and L z vs. Q r — (right-hand panel) space located inside a 
sphere of 2.5 kpc radius at 8 kpc from the galactic centre, after convolution with the Gaia observational errors. The black dots denote 
the contribution from the disc and bulge. The rest of the colours represent different satellites. 



4.2 Frequency Space 

Figure [9] shows the distribution of stellar particles inside the 
'Solar Neighbourhood' sphere in both frequency (left-hand 
panel) and L z vs. Q r — fl^ (right-hand panel) space, now 
after error convolution. Again, the different colours indicate 
stellar particles from different satellites, with the addition of 
the disc and bulge which are shown in black. Interestingly, in 
both spaces disc particles are part of a very well defined and 
quite isolated structure since, as expected for stars moving 
in the galactic plane on a circular orbit, they have the largest 
values of \L Z \ at a given il r — Q$. Therefore, it is possible to 
isolate and easily eliminate the contamination from the disc. 
On the other hand, bulge particles are sparsely distributed in 
both spaces and, although there are very few, they can not 
be simply isolated as in the case of the disc. Nonetheless, 
they do not define a clump that could be confused with 
merger debris. 

Both panels of Figure [9] show that, even after convolu- 
tion with the latest model of the Gaia measurement errors, 
significant substructure is apparent in L z vs. Q r — Q$ vs. Q$ 
space. Comparison to Figure [7] shows very good correspon- 
dence between clumps, although these are, as a consequence 
of the convolution, generally less well defined and contain 
less internal substructure. 

4-2.1 Estimating the time since accretion 

A derivation of the time since accretion is significantly more 
difficult when the Gaia measurement errors are taken into 
account. We exemplify this by focusing on the three satel- 
lites described in Section [3.2.21 The 'observed' distribution 
of stellar particles in frequency space are shown in grey in 
the top panels of Figure [TO] A direct comparison to Figure [8] 
clearly highlights what the effects of the Gaia errors are. As 
in Section [3.2.21 we perform our Fourier analysis and com- 
pute the normalized power spectra. These are shown with a 



black solid line in the bottom panels of Figure [POl Only for 
the spectrum shown in the middle panel, we can identify a 
single dominant peak. Therefore observational errors seem 
to be large enough to erase the signal associated to the time 
since accretion in the power spectrum for at least some of 
our satellite galaxies. 

We wish to obtain an order of magnitude estimate of 
how large the errors in velocities have to be to blend two 
adjacent streams in frequency space. Let us consider a satel- 
lite moving on a circular orbit accreted t a cc = 8 Gyr ago. 
The separation of two adjacent streams in frequency space 
is Af2 = 2ir/t acc w 0.78 (comparable to that found for 
the satellites in our simulations). Since V$ = RQ, we may 
deduce that the maximum error in the tangential velocity 
should be a V( , = RAQ. Therefore, for the streams found at 
R — 8 kpc from the Galactic centre, a Vtf> ^ 6 km s . This 
implies that, to be able to estimate the time since accretion 
from the power spectra, the relative parallax errors should 
be a^jvj < O.Ofl 

The set of stellar particles from each of our three satel- 
lites satisfying this condition are shown as red dots in the top 
panel of Figure ITUl Out of the 445 (left panel), 1264 (mid- 
dle panel) and 720 (right panel) stellar particles originally 
found inside our solar neighbourhood sphere, only 90, 239 
and 113 respectively, have remained after the error cut. The 
normalized power spectra obtained for these distributions 
are shown with dashed black lines in bottom panels of Fig- 
ure [10] Now, the largest amplitude peak is (once again and 
for all satellites) associated to the time since accretion. Note 
that, since each peak is located at the same wavenumber as 
in the analysis with no observational errors (grey curves in 



3 In this derivation we have assumed that relative errors in proper 
motion are of the same order of magnitude than those in the 
parallax. However, for the Gaia mission the former are expected 
to be in general an order of magnitude smaller. 
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Figure 10. Top panels: Distribution of stellar particles in frequency space for the three satellites shown in Fig. 8, after convolution 
with the Gaia observational errors. Grey dots show the distribution of all the stellar particles inside the 'Solar Neighbourhood' sphere, 
whereas the red dots are the subset with (r^/ro ^ 0.02 (see also zoom-in). Bottom panels: The black solid lines show the 1-D normalized 
power spectra obtained from the distribution of points shown in the top panels, while the black dashed lines corresponds to the subset 
with a^/zu ^ 0.02. The grey solid lines show the normalized power spectra obtained from the distribution of particles before convolution 
with the observational errors (cf. Figure [8}. Their largest amplitude peak is indicated whith a vertical line. 



this Figure), the estimated times since accretion for each 
satellite are exactly those obtained in Section \3. 2. 2 1 

Nevertheless, it is clear that the signal found in the 
power spectra could be determined better if a larger number 
of stars with a m /zu ^ 0.02 could be available. Such stars are 
expected to be present in the Gaia catalogue, but because 
of the limited numerical resolution of our experiments, they 
are not part of our Mock Gaia catalogue. Such stars would 
be fainter than Mv = 4.5, but closer than 2 kpc from the 
Sun. For example, according to the latest Gaia performance 
numbers, a dwarf star located at 1 kpc from the Sun, with 
an apparent magnitude of V ~ 16 should have a parallax 
measurement error a m /zu ^ 0.02. This apparent magnitude 
and distance translates into an absolute magnitude of My ~ 
6, i.e. it is 1.5 magnitudes fainter than limiting magnitude 
we have used in our simulations. If we use the number of 
particles found at distances closer than 1 kpc from the Sun, 
and if we take into account the luminosity functions of the 3 
satellites in Fig. 10, we can obtain an estimate of how many 
stars would be observable by Gaia with the desired parallax 
errors. In this way we find, respectively, 186, 465, and 276 
extra stars with 4.5 ^ Mv ^ 12, which will thus allows to 
estimate the time since accretion for each satellite. 



4.3 Identification of Satellites 

Although promising, this is unlikely to be the way we will 
proceed in the future with real Gaia observations. In or- 
der to obtain an estimate of the time since accretion of 
any given satellite, we first have to efficiently identify it. 
This can be performed by applying a suitable clustering 
technique . In this work we hav e used the Mean Shift al- 
gorithm ([F ukunag a fc Hostetlerl 1 19751 ; IComaniciu fc Meet] 
l2002i;lDerpanisll2005f ). The main idea behind mean shift is 
to treat the points in any iV-dimensional space as an em- 
pirical probability density function where dense regions in 
the space correspond to the local maxima of the underlying 
distribution. For each data point in the space, one performs 
a gradient ascent procedure on the local estimated density 
until convergence is reached. Furthermore, the data points 
associated with the same stationary point are considered 
members of the same cluster. 

We have applied the Mean Shift algorithm to the set 
of particles inside a sphere of 4 kpc radius, located at 
8 kpc from the galactic centre, projected into the space of 
E — L — L z . We chose this space because the satellite's inter- 
nal substructure (due to the presence of individual streams) 
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is less well defined and, therefore, is more suitable for a clus- 
tering search of global structures. Enclosed in this sphere we 
find a total of ~ 8 x f 4 stellar particles, coming from 26 dif- 
ferent satellites contributing each with at least 20 particles. 
In addition, the disc and the bulge contribute with ~ 20 000 
and ~ 2 000 particles, respectively. For this analysis we es- 
pecially chose to sample a larger volume of physical space 
so that each satellite is represented with a larger number of 
stellar particles. In this way, we can avoid overfragmenta- 
tion, which occurs when a clump is populated with a very 
small number of particles. 

Figure [TT] shows the different clusters of particles identi- 
fied with this method. We have found 17 groups that contain 
more than 20 particles. Out these 17 groups, only 15 have 
more than 50% of its particles associated to a single progen- 
itor. One of these groups corresponds to the disc while two 
other are double detections of two different satellites that 
were fragmented by the algorithm. Therefore, only 12 of 
these groups can be uniquely associated to a single satellite. 
This corresponds to ~ 50% of all the satellites contributing 
with stars to this volume. This is a s i milar recovery rate 
to that obtained bv lHelmi fc de Zeeuwl (|200d ), but now un- 
der a more realistic cosmological model and with the latest 
model for the Gaia measurement errors. 

When attempting to compute the time since accretion, 
we were successful only in four cases. These groups are in- 
dicated in Figure [11] with black open circles. Two of them 
correspond to the satellites labeled number 1 and 2 in previ- 
ous analysis. In general, we find that the remaining identified 
groups lack a significant number of 'stars' with the required 
relative parallax error, a^/vj ^ 0.02 (i.e., typically ^ 50). 
However, as explained before, in this simulation our lim- 
ited numerical resolution led us to consider only stars with 
Mv ^ 4.5. After estimating the number of stellar particles 
within 1 kpc from the Sun that may be observed by Gaia 
with the required relative parallax errors (as explained in 
Section I4.2.ip . we find that at least two additional satel- 
lites, among the 12 previously isolated, should have at least 
200 stellar particles available to compute a reliable power 
spectrum. 



5 SUMMARY AND CONCLUSIONS 

We have studied the characteristics of merger debris in the 
Solar neighbourhood as may be observed by ESA's Gaia mis- 
sion in the near future. We have run a suite of JV-body sim- 
ulations of the formation of the Milky Way stellar halo set 
up to match, at the present time, its known properties such 
as the velocity ellipsoid, density profile and total luminosity. 
The simulations follow the accretion of a set of 42 satellite 
galaxies onto a semi-cosmological time dependent Galactic 
potential. These satellites are evolved for 10 Gyr, and we use 
the final positions and velocities of the constituent particles 
to generate a Mock Gaia catalogue. 

Using synthetic CMDs, we have resampled the satel- 
lite's particles to represent stars down to My « 4.5. 
This absolute magnitude corresponds to an apparent mag- 
nitude V ~ 16 at 2.5 kpc, which is comparable to the 
Gaia magnitude limit for which full phase-space informa- 
tion will be available. Our Mock catalogue also includes 
a Galactic background population of stars represented by 



x 10 




-3000 -2000 -1 000 



1000 2000 3000 



L [kpc km s ] 

Figure 11. Distribution of stellar particles inside a sphere of 
4 kpc radius at 8 kpc from the galactic centre in E vs. L z space 
as would be observed by Gaia. The different colours show the 
groups identified by the Mean Shift algorithm. Black open circles 
denote those for which the time since accretion was successfully 
derived. 



a Monte Carlo model of the Galacti c disc and bulge, as 
in IfBrown. Velasquez fc Aguilarl l|2005l h At 8 kpc from the 
Galactic centre, stars from the disc largely outnumber those 
of the stellar halo, however it is possible to reduce their 
impact by applying a simple cut on metallicity. We have es- 
timated, using t he latest determinations of the MDF of th e 
Galactic disc(s) (|Nordstrom et al.l 1200-1 llvezic et alj|2008l ). 
that only 10 000 stars out of the estimated 4.1 x 10 7 disc stars 
brighter than V = 17 in the Solar neighbourhood should 
have [Fe/H] ^ —1.1. A smaller number of bulge stars (~ 800) 
with [Fe/H] ^ —1.1 is expected to contaminate our stellar 
halo sample, down to the faintest My. This fraction repre- 
sents only 2.6 per cent of the whole Mock Gaia stellar halo 
catalogue, and therefore does not constitute an important 
obstacle to our ability to characterize this component. 

Finally, we have convolved the positions and veloc- 
ities of all 'observable stars' in our Solar Neighbour- 
hood sphere with the latest models of the Gaia obser- 
vational errors, according to performances given by ESA 
(http:/ /www.rssd.esa.int/gaia). 

The analysis presented here confirms previous re- 
sults, namely that satellites can be identified as coher- 
ent clumps in phase-space, e.g. in the E-L z projection 
fseelHelmi fc de Zeeuwl hood : iKnebe et al1l2005l ; iFont et all 
2006) . We find that a clustering algorithm such as Mean- 
Shift (|Fukunaga fc Hostetlerlll975l ; IComaniciu fc Meej2002l : 



|Perpanisl l2005) is able to recover roughly 50% of all satel- 
lites contributing stellar particles to the Solar neighbour- 
hood sphere. 

We have also demonstrated that even after accounting 
for the Galactic background contamination and the expected 
measurement errors for the Gaia mission, the space of or- 
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bital frequencies is also very rich in substructure. In this 
space streams from a given satellite define a regular pattern 
whose characteristic scale is determined by the satellite's 
time since accretion. We find that reasonable estimates of 
the time since accretion may be provided when the num- 
ber of stars with accurate parallaxes {a-^/zu ^ 0.02) from 
a given satellite is at least ~ 100. This was possible in our 
simulations for 4 cases. 

Together with Gaia, ground based follow-up campaigns 
are currently being planned by European astronomical com- 
munity. Their purpose is to complement a future Gaia cat- 
alogue with information that either will not be obtained 
or for which the accuracy will be very low. As an exam- 
ple, high resolution multi object spectroscopy in combina- 
tion with ground based photometry could not only push 
the limiting magnitude of the phase-space catalogue down 
to V ~ 20, but could also provide detailed chemical abun- 
dance patterns of individual stars. The former will require 
the development of accurate and precise photometric dis- 
tances indicators which can be calibrated using the large 
number of stars near the Sun with very precise trigono- 
metric parallaxes from Gaia. The identification of satellites 
in, e.g., E — L z space could be considerably improved by 
having such an extended full phase-space catalogue. Firstly, 
very faint satellites could now be observed and, secondly, 
the much larger sample of stars could allow us to apply 
a clustering algorithm in very small volumes around the 
Sun with enough resolution to avoid large overfragmenta- 
tion. The later is important because by extending the vol- 
ume probed we expect to reduce the overlap between satel- 
lites in E — L z space. In addition, chemical tagging will be 
very important to char acterize the satellite's star formatio n 
and chemical histories (|Freeman fe Bland- Hawthorol |2002|) . 
Note, however, that in our simulations we find that stellar 
particles in a given stream do not originate from a localized 
region in physical space (such as a single molecular cloud). 
Therefore even individual streams are likely to reflect the 
full MDF present in the object at the time it was accreted. 

Although our satellites were evolved in a cosmologically 
motivated time dependent host potential, our simulations do 
not contain the fully hierarchical and often chaotic build-up 
of structure characteristic of the A cold dark matter model. 
The violent variation of the host potential during merger 
events, the c haotic behavior induced b y a triaxial dark mat- 
ter halo (e.g. IVogelsberger et al.ll200Sl), and the orbital evo- 
lution due to baryonic condensation (|Valluri et alj20ich , are 
potentially important effects which we have neglected and 
should be taken into account in future work. To assess the 
impact of some of these effects on the distribution of de- 
bris in the Solar Neighbourhood we are currently analysing 
a fully cosmological high resolution simulation of the for- 
mation of the Galactic stel lar halo based on the Aquarius 
project (|Cooper et alj|2010h . 
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